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A new Monte Carlo method for long-range interacting systems is presented. This method 
involves eliminating interactions stochastically with the detailed balance condition satisfied. 
When pairwise interactions Vij of an A^-particle system decrease with the distance as r~j°', 
computational time per Monte Carlo step is 0{N) for a > d and 0{N^~°'^'^) for a < d, where 
d is the spatial dimension. We apply the method to a two-dimensional magnetic dipolar 
system. The method enables us to treat a huge system of 256^ spins within a reasonable 
computational time, and reproduces a circular order originating from long-range dipolar 
interactions. 
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Numerical simulations of long-range interacting systems are quite difficult because we have 
to take a large number of interactions into consideration. If one carries out a naive Monte 
Carlo (MC) simulation of an A^-particle system with pairwise interactions, computational time 
per MC step tMC is proportional to A^^. Due to this rapid increase in computational time, 
accessible sizes for numerical simulations are restricted. One might think that this problem 
can be resolved by truncating interactions beyond a certain cutoff distance. However, such 
truncations often bring significant errors in various observables.^^^^ Concerning long-range 
interacting Ising ferromagnetism, a cluster algorithm that drastically improves computational 
efficiency without any approximation has been proposed.^) However, this method cannot be 
used for other systems. To overcome this difficulty, some approximate methods have been 
proposed until now.®"^^^ Some of them can treat more than one million particles within a 
reasonable computational time with high accuracy. Furthermore, these methods are appli- 
cable to general long-range interacting systems. Nevertheless, these methods include some 
approximations more or less. 

In this paper, we present a new MC method for general long-range interacting systems. In 
contrast to other methods, the present method is exact in the sense that it strictly satisfies the 
detailed balance condition. In this method, we stochastically switch long-range interactions 
Vij to either zero or a pseudointeraction Vij by use of the stochastic potential switching 
algorithm. Then the system is mapped on that with only Vij. The potential switching is 
performed every several steps. Since most of the distant (and weak) interactions are eliminated 
(Vij 0), tuc is significantly reduced. We refer to the present method as the stochastic cutoff 
(SCO) method. Of course, if one naively switches Vij, it costs computational time tswitch of 
order N'^. We develop an efficient method for the potential switching. In lattice systems, 
it reduces tgwitch to be comparable to tuc- We apply the SCO method to a two-dimensional 
magnetic dipolar system. By comparing our data with the previous ones,^^^ we confirm that the 
SCO method gives correct results with modest computational time. We discuss the properties 
of the SCO method in comparison with other methods for long-range interacting systems. 

The organization of the paper is as follows. In § 2, we describe the SCO method. In § 3, 
we show the results obtained by applying the SCO method to a two-dimensional magnetic 
dipolar system. In § 4, we compare the SCO method with other methods. Section 5 is devoted 
to conclusions. 

2. Stochastic Cutoff" (SCO) Method 

Before explaining the SCO method, we briefly survey the stochastic potential switching 
algorithm. We hereafter consider a lattice system with pairwise long-range interactions 
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described by the Hamiltonian 

H = Y,Vi,{S^,Sj), (1) 

where Si is a variable associated with the i-th element of the system. In this algorithm, Vij 
is stochastically switched to either Vij or Vij with a probability of Pij or 1 — Pjj, respectively. 
The probability Pij is 

Pij{S,, S,) = eMP{/Wij{Si, Sj) - AV^j)], (2) 

where (3 is the inverse temperature, AVij{Si, Sj) = Vij{Si, Sj) — Vij{Si, Sj), and AV*j is a 
constant equal to (or greater than) the maximum value of AVij{Si, Sj) over all Si and Sj. We 
can choose the potential Vij arbitrarily. On the other hand, using Pij{Si, Sj), the potential 
Vij is given as 

VjiSi, S,) = VjiSi, Sj) - log[l - P^J{S„ Sj)]. (3) 

The arbitrariness of Vij can be utilized to reduce either the complexity of potential or cost to 
calculate it. With this potential switching process, the algorithm proceeds as follows: 

(A) Potentials Vij are switched to either Vij or Vij with the probability of Pij or 1 — Pij, 
respectively. 

(B) A standard MC simulation is performed with the switched Hamiltonian 7i' expressed as 

where runs over all the potentials switched to V and runs over those switched to 
V. The potential is fixed during the simulation. 

(C) Return to (A). 

It is shown that this MC procedure strictly satisfies the detailed balance condition with respect 
to the original Hamiltonian of eq. (1). 

We give two remarks. Firstly, we can choose the period of the simulation in step (B) 
arbitrarily because the detailed balance condition is satisfied regardless of the period. Secondly, 
we do not need to switch all the potentials. For example, when the original Hamiltonian is 
^ — ^ij + ^ij^ we can switch only {Vij} with {Uij} unswitched. Such partial switching 
is realized by setting Uij = Uij. 

We now describe the SCO method. The basic idea is quite simple. We just set Vij = to 
reduce computational time. In the following, we see that time tuc in step (B) is significantly 
reduced since most potentials associated with distant interactions are switched to Vij = in 
step (A). Now let us assume that interactions Vij decrease as Dfij"' > 0), where Z) is a 
constant that represents the strength of interactions and rij is the distance between sites i 
and j. Systems with dipolar interactions correspond to the case a = 3. When we update an 
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element from to sf^™^ in a MC simulation, we need to calculate the energy difference 

" El [^^'^ (^S""^' ^'^) - {st'^^S,)] , (5) 
where the sum is taken over all the sites k for which Vik is switched to Vik- This means that 
the time required to update Si is proportional to the number of terms, i.e.. Mi = When 
Vik is large, the probability that Vik is switched to Vik is approximated as 



1-P, 



ik 



D(3rr- (ra. » (I?/?)'/") , (6) 



where we have assumed that both AVik and AT^^ are of order Dr^^. Therefore, we can 
roughly evaluate Ni as 



tlAC oc NMi ^ < 



{a > d), 

D(3\og{L) {a = d), (7) 
D(3L'^~'' {a<d), 

where d is the spatial dimension of the system, L is the linear size of the lattice, and the 
lattice constant is assumed to be one. Therefore, computational time per one MC step ^mc is 
estimated as 

D(3N {a > d), 
D/3Nlog{N) {a = d), (8) 

where N{= L'^) is the total number of elements and the difference between log(L) and log(A^) 
is ignored. Obviously, t^c in the present method is much smaller than that in a naive MC 
method where tuc ~ 

We next consider the potential switching process. If one switches Vij one by one, the 
potential switching time tswitch is of order A^^. To reduce tswitch) we develop the following 
method. We first introduce a set of pairs {Cpairl'")} for which either Vij or Vji is r, where Vij 
is the vector spanning from sites i to j. Figure 1 shows an example of such sets. Pairs are 
labeled sequentially. The point is that, in most cases, the probabilities of switching to V for 
{Cpaiiir)} have some upper limit PmaxC'")- As we will show later, such maximum probability 
indeed exists in dipolar systems. Using Pma.x{f), we can switch the potentials in {Cpair(T')} in 
the following way: 

a) Using pmax(^)) we choose candidates that are switched to Vij. The potentials that are not 
chosen candidate are switched to 0. 



b) Switch each of the candidates to Vij with the probability [1 — Pij{Si, Sj)]/pinsLx{r). Oth- 
erwise, Vij is switched to 0. 
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Of course, if candidates are chosen one by one, it is very timeconsuming. We therefore choose 
them in the fohowing way. We hereafter denote the potential of the k-th pair in {Cpair(*')} 
by V^'^l Since the probabihty of being chosen as a candidate is the same in {Cpair the 
probabihty that V^^^ is chosen as a candidate after n — 1 successive failures is given by a 
geometric distribution g{pma.x{f),n), where 

g{p,n) = {l-pr-'p (n>l). (9) 

An integer random variate n that obeys g{p, n) can be easily generated as 

„=rj2i(!i^i, (10) 

log(l - p) 

where \x\ is the smallest integer that is greater than or equal to x, and r is a continuous 
random variate with an uniform distribution of range < r < 1. We can pick up only 
candidates by means of fl'(pmax('")) '^)- For example, the generation of two random variates n\ 
and n2 means that there are only two candidates y('^i) and y("i+"2) among ni+n2 potentials. 
Using this idea, we have implemented the potential switching in {CpairC^")} as follows: 

1) Set Tig to be 0, where Ug is the number of potentials that have already been switched. 

2) Generate an integer n from the distribution g{praaxi'r'),n) using eq. (10). If n = 1, go to 
step 4). Otherwise, go to the next step. 

3) Switch the n - 1 potentials (F^+i), . . . , -^K+n-i)) to 0. 

4) Switch to t/K+") with the probability [1 - p(''=+")(S'i, Sj)]/pm!,Ar)- Otherwise, 
switch it to 0. 

5) Finish the potential switching procedure if ris + n is greater than (or equal to) the number 
of elements of {Cpair('')}- Otherwise, replace with Tig + n and return to step 2). 

Switching of all the potentials is completed by carrying out this procedure for all {Cpair ('")}• 
Now let us evaluate ^switch- The potential switching process clearly requires time propor- 
tional to the number of potentials chosen as a candidate. To estimate it, we focus on the A^ — 1 
potentials associated with a certain site i and estimate the number of candidates Mj^^^ among 
them. The total number of candidates is AA/'/'^^ Since both Pma.x{f) and 1 — Pj^ are of order 
D(3r~°' when r ^ {Dj3Y^°^ , the order estimation of M^"'^ and that of Mi given by eq. (7) are 
the same. We therefore obtain 

D^N {a > d), 
tswitch~< D(3NlogiN) {a = d), (11) 

From eqs. (8) and (11), we find that tgwitch is indeed comparable to tMC- 
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Now we apply the SCO method to a two-dimensional magnetic dipolar system on an L x L 
square lattice with open boundaries. The Hamiltonian of the system is described as 

n = -jy^ Si -Si 



i<j 



(12) 

where Si is a classical Heisenberg spin of \Si\ = 1, {ij) runs over all the nearest-neighboring 
pairs, Vij is the vector spanned from sites i to j in the unit of the lattice constant a, and 
fij = \ The first term describes short-range ferromagnetic exchange interactions and the 
second term describes long-range dipolar interactions, where J(> 0) is an exchange constant 
and D = {gii^S)/a?. Hereafter, we regard D as a parameter and consider the case that 
D/ J = 0.1. We choose this model as a benchmark of the SCO method because the properties of 
the model have been investigated extensively in previous work.^^^ In particular, it is established 
that the model undergoes a phase transition from the paramagnetic state to a circularly 
ordered state at ~ 0.88J as a consequence of the cooperation of exchange and dipolar 
interactions. 

Before showing the results, we explain the details of our simulation. We applied the SCO 
method only for dipolar interactions. The potential difference AV^j is given as 



mj = = D 



(13) 



It has the minimum value —2D/rfj and the maximum value +2D/rfj when Si and Sj are 
parallel along Vij and antiparallel, respectively. Therefore, we obtain 

= max^^ 5^.{l - exp[P{AVi,{S„ Sj) - AV*j)]] 

= 1 -exp[-4D/3/r2], (14) 

where we have set AV*j = 2D/r^. The system was gradually cooled from an initial temperature 
T = 1.65J to 0.05J in steps of AT = 0.05J. The initial temperature was set to be well above 
the critical temperature. The system was kept at each temperature for 100, 000 MC steps, and 
potentials were switched for every 10 MC steps. The first 50, 000 MC steps are for equilibration 
and the following 50,000 MC steps are for measurement. Therefore, the total MC steps for 
one run is 3, 300, 000. The computational time per run for L = 256, i.e, the maximum size we 
examined, was less than three days when using a personal computer with a Core2Duo 2.40 
GHz processor. We conducted simulations for 10 different runs with different initial conditions 
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and random sequences. The energy was calculated for every 10 MC steps with 0{N\ogN) 
computational time by utilizing the discrete Fourier convolution theorem and the fast Fourier 
transformation algorithm. For a detailed description of how the discrete Fourier convolution 
theorem is used for a system with open boundary conditions, we refer the reader to ref. 14. 

Now let us see the results of our simulations. Figure 2 shows the temperature dependences 
of the specific heat for different sizes. The peaks are located around Tc ~ 0.88J. In Fig. 3, 
we show the spin structure observed at T = 0.05 J. We clearly see a circular order that comes 
from long-range dipolar interactions. For a quantitative measurement of the circular order, 
we observed the absolute value of the circular component defined by 



(15) 



where (• • • ) denotes the thermal average and Vc is a vector describing the center of the lattice. 
Figure 4 shows the result. We see that the circular order rapidly grows around the critical 
temperature. We also performed a naive MC simulation for L = 48 to confirm that the SCO 
method reproduces correct results. The crosses in Figs. 2 and 4 show the results. 

Now let us examine the efficiency of the SCO method. Figure 5 shows the size dependences 
of the average computational time per MC step tav for both the SCO method and the naive 
MC method. The average time tav of the SCO method is given as 

^av = ^MC + ^Q^switch + Yg^energy, (16) 

where tenergy is the computational time per one energy measurement. Recall that potential 
switching and energy measurement are performed for every 10 MC steps. We see that tav oc N 
in the SCO method, which is strongly in contrast with tav oc N"^ in the naive MC method. 

This huge reduction of the computational time in the SCO method comes from the reduc- 
tion of interactions. We observed the average number n of potentials per site that survive as 
Vij. In Fig. 6, we show the temperature dependences of n for different sizes. It is impressive 
that n at each temperature converges to a certain value as the size increases. Although n 
increases with decreasing temperature, n ~ 22.5 even at the lowest temperature. 

In Fig. 7, we show the data of the specific heat measured in several simulations with 
different potential switching periods N^v,- Since the system is kept at each temperature for 
100,000 MC steps, the number of potential switchings at each temperature is 100, OOO/A'^s^- 
We see that reliable results are obtained when N^^ is 1,000 or less. Figure 8 shows the size 
dependence of the ratio tswitch/^MC- The ratio slightly depends on the size, and it is 1.65 at 
most. This means that we can even switch potentials at every MC step with a reasonable cost. 

Lastly, we compare relaxation speeds between the SCO method and the naive MC method. 
Figure 9 shows the time evolution of for L = 48 when the system is kept at T = 0.4J. As 
shown in the inset, the relaxation speed of the SCO method is about 1.4 times slower than 
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that of the naive MC method. We have also performed a similar measurement at T = 0.7 J, 
and found that the ratio is about 1.2. It is clear from Fig. 5 that the relaxation speed of the 
SCO method is much higher than that of the naive MC method if they are compared in terms 
of the computational time. 

4. Comparison of the SCO Method v^ith Other Methods 

In this section, we discuss the properties of the SCO method in comparison with those 
of other methods for long-range interacting systems. As we have emphasized, the primary 
merit of the SCO method is that it involves no approximation. The performance of this 
method strongly depends on the conditions of simulations in terms of parameters such as the 
temperature T, the spatial dimension d, the decay exponent of potentials a, and the strength 
of potentials D, as illustrated by eqs. (8) and (11). Since tMC and tgwitch are proportional to 
Dp, the SCO method is particularly efficient for systems with strong short-range interactions 
and weak long-range interactions. The reason is as follows: Since the system is dominated by 
short-range interactions, we expect /ebTc ~ -^sr, where -DsR is the strength of short-range 
interactions. This means that tMC and tswitch are very small around the critical temperature 
because 

D(3c ~ D/DsR < 1. (17) 

In this sense, the SCO method is really suitable for magnetic dipolar systems, as we have 
demonstrated. In magnetic dipolar systems, the ratio J/D is usually on the order of hundreds 
or even thousands. On the other hand, when a < d, the computational time in the SCO method 
is 0{N'^~'^^'^). The SCO method is less efficient in such a case because the computational time 
is of order N or N log N in most of the other methods. 

Lastly, we emphasize that the SCO method can be made applicable to off-lattice systems 
by developing some efficient method for potential switching. 

5. Conclusion Remarks 

In the present work, we have proposed a new MC method based on the stochastic potential 
switching algorithm. ^^'^^^ To our knowledge, this is the first method for general long-range 
interacting systems that greatly reduces computational time without any approximation. This 
method is applicable to any lattice system with long-range interactions. The efficiency of the 
SCO method has been demonstrated by applying it to a two-dimensional magnetic dipolar 
system. We have also discussed the properties of the SCO method in comparison with those 
of other methods for long-range interacting systems. 
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Figure 1 : Set of pairs {Cpair{r)} (r = (1, 1)) and their labels in a 5 x 5 square lattice. The set 
consists of 16 elements. 

Figure 2 : (Color online) Temperature dependences of the specific heat C for different sizes. 
The average is taken over 10 different runs. 

Figure 3 : Snapshot of the spin structure at T = 0.05J on a 48 x 48 square lattice obtained 
by the SCO method. 

Figure 4 : (Color online) Temperature dependences of the circular component of the magne- 
tization for different sizes. The average is taken over 10 different runs. 

Figure 5 : (Color online) Average computational time per MC step tav is plotted as a function 
oi N = LP' for the SCO method (full squares) and a naive MC method (full circles). The solid 
line and dashed line are proportional to x and x^, respectively. 

Figure 6 : (Color online) Temperature dependences of the average number n of potentials 
iyij 7^ 0) per site. The average is taken over 10 different runs. 

Figure 7 : (Color online) Temperature dependences of the specific heat on a 48 x 48 square 
lattice with different N^-w values, where A'g^ is the number of MC steps for every which po- 
tentials are switched. The data measured in a single run are shown. 

Figure 8 : Ratio tswitch/^MC is plotted as a function of = L^. 

Figure 9 : (Color online) Time evolution of the circular component of the magnetization 
for L = 48 when the system is kept at T = 0.4J. The average is taken over 1,000 different 
runs. The data of the SCO method are denoted by crosses (below) and those of the naive MC 
method are denoted by open squares (above). In the inset, of the SCO method and that 
of the naive MC method are plotted as a function oi = t/lA and t, respectively. 
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